module count
    implicit none
    save
    integer(4) :: jsn,jcn
end module count

program xelf
!     test driver for "elf"
    use complete_elliptic
    implicit none
    !real*8 PI,PIHALF
    !parameter (PI=3.1415926535897932384626433d0)
    !parameter (PIHALF=1.5707963267948966192313216916398d0)
    real*8 dmc,dphi,mc,mm,phi,phic,xf,elf
    integer jend,iend,j,i
    !integer,parameter :: dp=8

    jend=6
    iend=5
    dmc=1._dp/dble(jend-1)
    dphi=PIHALF/dble(iend)
    write(*,'(1x,2a10,a25,a15)')	'm','phi/PI','elf','relf'
    do j=1,jend
        write(*,'(1x)')
        mc=real(j-1,dp)*dmc
        if(mc.le.0._dp) mc=1.21e-32_dp
        mm=1._dp-mc
        do i=0,iend
            phi=dphi*real(i,dp)
            phic=dphi*real(iend-i,dp)
            xf=elf(phi,phic,mc)
            write(*,'(1x,0p2f10.5,0p1f25.15)') mm,phi/PI,xf
        enddo
    enddo
    END
!---------------------------------------------------------------------------
    function elf( phi , phic , mc )
!
!	Double precision incomplete elliptic integral of the first kind
!
!     Reference: T. Fukushima, (2010) Numer. Math. 116, 687-719
!        "Fast Computation of Incomplete Elliptic Integral of First Kind
!         by Half Argument Transformation"		

!     Author: T. Fukushima Toshio.Fukushima@nao.ac.jp

!     Used subprograms: asn, acn, elk, serf (called from asn)

!     Inputs: phi  = argument                0 <= phi  <= PI/2
!             phic = complementar argument   0 <= phic <= PI/2
!             mc   = complementary parameter 0 <= mc   <= 1

!     Output: elf

!     CAUTION: phi and phic must satisfy condition, phi + phic = PI/2
    use complete_elliptic
    implicit none
    !integer,parameter :: dp=8

    real(kind=dp),intent(in) :: phi,phic,mc
    real(kind=dp) asn,acn,m,c,yc,d2,v,elf,CompleteEllipticFirstKind

    m=1._dp-mc
    
    if ( phi < 1.25_dp ) then
        elf = asn( sin(phi) , m )
    else
        c = sin(phic)
        yc = c*c
        d2 = mc+m*yc
        if ( yc < 0.9_dp*d2 ) then
            elf = CompleteEllipticFirstKind(mc) - asn( c/sqrt(d2) , m )
        else
            v = mc*(1._dp - yc)
            if ( v < yc*d2 ) then
                elf = acn( c , mc )
            else
                elf = CompleteEllipticFirstKind(mc) - asn( sqrt(v/d2) , m )
            endif
        endif
    endif
    return
    end function elf


    function acn(c0,mc)
    use count
    implicit none
    integer,parameter :: dp=8

        real(kind=dp),intent(in)        :: c0,mc
        real(kind=dp)                   :: m,c,yc,y,s,f,d,asn,acn
    
        integer   ::i    !,jsn,jcn
        !common /count/ jsn,jcn

        m  = 1._dp - mc
        c  = c0
        yc = c*c
        
        if ( yc > 0.5_dp ) then
            y = 1._dp - yc
            s = sqrt(y)
            acn = asn(s,m)
            jcn = 0
            return
        endif
        
        f = 1._dp
        
        do i=1,20
            d = sqrt( mc + m*yc )
            yc = (c + d)/(1._dp + d)
            f = f*2._dp
            
            if( yc > 0.5_dp ) then
                y = 1._dp - yc
                s = sqrt(y)
                jcn = i
                goto 1
            endif
            
        c = sqrt(yc)
        enddo
        
        write(*,*) "(acn) too many iterations: c0,mc=",c0,mc
    1   continue
        acn = f*asn(s,m)
        return
    end function acn

    function asn(s0,m)
    use count
    implicit none
    integer,parameter :: dp=8

        real(kind=dp),intent(in)        :: s0,m
        real(kind=dp)                   :: del,s,f,y,serf,asn
        integer                             :: j                        !jsn,jcn , j
        !common /count/ jsn,jcn

        del = 0.04094_dp - 0.00652_dp*m
        s = s0
        y = s*s
        if ( y < del ) then
            asn = s*serf(y,m)
            jsn = 0
        return
        endif
        f = 1._dp
        do j=1,20
            y = y/( (1._dp+sqrt(1._dp-y))*(1._dp+sqrt(1._dp-m*y)) )
            f = f*2._dp
            if ( y < del ) then
                s = sqrt(y)
                jsn = j
                goto 1
            endif
        enddo
        write(*,*) "(asn) too many iterations: s0,m=",s0,m
    1   continue
        asn = f*s*serf(y,m)
        return
    end function asn

    function serf(y,m)
    implicit none
    integer,parameter :: dp=8

        real(kind=dp),intent(in) ::y,m

        real(kind=dp)                      :: F1,F2,F3,F4,F5,F6,F7,F8,F9,serf
        real(kind=dp),parameter     ::F10=1._dp/6._dp
        real(kind=dp),parameter     ::F20=3._dp/40._dp
        real(kind=dp),parameter     ::F21=2._dp/40._dp
        real(kind=dp),parameter     ::F30=5._dp/112._dp
        real(kind=dp),parameter     ::F31=3._dp/112._dp
        real(kind=dp),parameter     ::F40=35._dp/1152._dp
        real(kind=dp),parameter     ::F41=20._dp/1152._dp
        real(kind=dp),parameter     ::F42=18._dp/1152._dp
        real(kind=dp),parameter     ::F50=63._dp/2816._dp
        real(kind=dp),parameter     ::F51=35._dp/2816._dp
        real(kind=dp),parameter     ::F52=30._dp/2816._dp
        real(kind=dp),parameter     ::F60=231._dp/13312._dp
        real(kind=dp),parameter     ::F61=126._dp/13312._dp
        real(kind=dp),parameter     ::F62=105._dp/13312._dp
        real(kind=dp),parameter     ::F63=100._dp/13312._dp
        real(kind=dp),parameter     ::F70=429._dp/30720._dp
        real(kind=dp),parameter     ::F71=231._dp/30720._dp
        real(kind=dp),parameter     ::       F72=189._dp/30720._dp
        real(kind=dp),parameter ::       F73=175._dp/30720._dp
        real(kind=dp),parameter ::      F80=6435._dp/557056._dp
        real(kind=dp),parameter ::      F81=3432._dp/557056._dp
        real(kind=dp),parameter ::      F82=2722._dp/557056._dp
        real(kind=dp),parameter ::      F83=2520._dp/557056._dp
        real(kind=dp),parameter ::      F84=2450._dp/557056._dp
        real(kind=dp),parameter ::    F90=12155._dp/1245184._dp
        real(kind=dp),parameter ::      F91=6435._dp/1245184._dp
        real(kind=dp),parameter ::      F92=5148._dp/1245184._dp
        real(kind=dp),parameter ::      F93=4620._dp/1245184._dp
        real(kind=dp),parameter ::      F94=4410._dp/1245184._dp

        F1=F10+m*F10
    
        F2=F20+m*(F21+m*F20)
    
        F3=F30+m*(F31+m*(F31+m*F30))
    
        F4=F40+m*(F41+m*(F42+m*(F41+m*F40)))
    
        F5=F50+m*(F51+m*(F52+m*(F52+m*(F51+m*F50))))
    
        F6=F60+m*(F61+m*(F62+m*(F63+m*(F62+m*(F61+m*F60)))))
    
        F7=F70+m*(F71+m*(F72+m*(F73+m*(F73+m*(F72+m*(F71+m*F70))))))
    
        F8=F80+m*(F81+m*(F82+m*(F83+m*(F84+m*(F83+m*(F82+m*(F81+m*F80)))))))
    
        F9=F90+m*(F91+m*(F92+m*(F93+m*(F94+m*(F94+m*(F93+m*(F92+m*(F91+m*F90))))))))
    
        serf=1._dp+y*(F1+y*(F2+y*(F3+y*(F4+y*(F5+y*(F6+y*(F7+y*(F8+y*F9))))))))
        return
    end function serf
